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Different numerical methods can be applied to the analysis of the flow through the Stirling engine regen¬ 
erator. One growing approach is to model the regenerator as porous medium to simulate and design the 
full Stirling engine in three-dimensional (3-D) manner. In general, the friction resistance coefficients and 
heat transfer coefficient are experimentally obtained to describe the flow and thermal non-equilibrium 
through a porous medium. A finite volume method (FVM) based non-thermal equilibrium porous media 
modelling approach characterizing the fluid flow and heat transfer in a representative small detailed flow 
domain of the woven wire matrix is proposed here to obtain the porous media coefficients without 
further requirement of experimental studies. The results are considered to be equivalent to those 
obtained from the detailed woven wire matrix for the pressure drop and heat transfer. Once the equiv¬ 
alence between the models is verified, this approach is extended to model oscillating regeneration cycles 
through a full size regenerator porous media for two different woven wire matrix configurations of 
stacked and wound types. The results suggest that the numerical modelling approach proposed here 
can be applied with confidence to model the regenerator as a porous media in the multi-dimensional 
(multi-D) simulations of Stirling engines. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

The Stirling engine is an attractive alternative machine for 
micro-cogeneration mainly for the achievable high efficiencies 
and the possibility to use different heat sources. Recent studies 
on Stirling engines are motivated by the need to reduce the energy 
consumption, to improve energy efficiency and to use clean energy 
technologies. Despite of having many theoretical advantages, Stir¬ 
ling engines still face some challenges in terms of efficiency and 
the design process. 

The recent Stirling engines studies mainly focus on the multi¬ 
dimensional (multi-D) analysis [1-7], which is based on full Stir¬ 
ling engine modelling, and can reproduce useful Stirling engine 
performance results in a time-frame short enough to impact design 
decisions [3 . In the multi-D analysis, the regenerator is a difficult 
component to be modelled because any numerical inaccuracies 
will influence the full scale Stirling simulation. In the modelling 
of a regenerator, an important factor is the internal geometry of 
the matrix and most of the regenerator models do not assume a 
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precise geometrical shape for the elements of the regenerator. In 
the Stirling engine’s multi-D analysis the regenerator is usually 
modelled as a macro-scale porous medium, the porous media 
model requires the input of a friction factor and heat transfer coef¬ 
ficient which mostly are empirically obtained. 

There are several works on modelling of Stirling regenerator as a 
porous media [8-16 . Cleveland State University 8] worked on CFD 
to develop a multi-D computer model of a portion of the regenera¬ 
tor matrix. Park et al. [9] indicated that screen-laminate matrices 
can be modelled as a porous media using a local thermal non-equi¬ 
librium model. Kim [10] studied the flow friction associated with 
laminar pulsating flows through porous media mainly for under¬ 
standing the Stirling regenerators and pulse tube cryocoolers. Tao 
et al. [11] investigated the fluid flow and heat transfer perfor¬ 
mances of mesh regenerators under different mesh geometric 
parameters and material properties based on an anisotropic porous 
media. Landrum et al. [12] modelled porous media hydrodynamic 
parameters adjusting iteratively to match the model predictions 
to the experimental results. Conrad et al. [13] determined the 
hydrodynamic parameters of stacked discs matrices using a 2-D 
CFD assisted methodology, whereby the hydrodynamic resistance 
parameters for these fillers are specified when they are modelled 
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Nomenclature 



Afs 

interfacial area density (m 2 ) 

Si 

source term (Pa/m) 

a\ , a2 

friction factor correlation equation constants (-) 

c h 

Jf.s 

fluid/Solid enthalpy source (W/m 3 ) 

Bf 

body forces (Pa/m) 

S 

cell size (m) 

c 2 

inertial resistance factor (-) 

Tf 

working gas temperature (K) 

C E 

Ergun coefficient (-) 

T s 

solid medium temperature (I<) 

Cf 

friction factor (-) 

t 

time (s) 

dh 

regenerator matrix hydraulic diameter (m) 

u 

inlet velocity magnitude in x direction (m/s) 

d w 

wire diameter (m) 

u 0 

amplitude of oscillating inlet velocity (m/s) 

Bf 

total fluid energy (J/kg) 

Umax 

matrix velocity magnitude in x direction (m/s) 

E s 

total solid energy (J/kg) 

V 

overall velocity vector (m/s) 

f 

frequency (Hz) 

a 

intrinsic permeability (m 2 ) 

h 

heat transfer coefficient (W/m 2 K) 

p 

Forchheimer coefficient (1/m) 

hf s 

heat transfer coefficient for fluid/solid interface (W/ 

P 

fluid dynamic viscosity (kg/m s) 


m 2 I<) 

Pf 

fluid density (kg/m 3 ) 

Bf 

working gas thermal conductivity (W/m I<) 

Ps 

solid density (kg/m 3 ) 

ks 

solid thermal conductivity (W/m 2 K) 

n v 

regenerator matrix volumetric porosity (-) 

L 

length of regenerator matrix (m) 

V 

regenerator matrix specific heat transfer area (1/m) 

Nu 

Nusselt number (-) 

y 

porosity of the media (-) 

P 

static pressure (Pa) 

T 

non-dimensional time (-) 

Vp 

pressure gradient (Pa/m) 

f 

stress tensor (Pa) 

Re 

Reynolds number (-) 




as anisotropic porous media. Nair and Krishnakumar [14] used a 
2-D local thermal equilibrium porous media model to investigate 
the heat transfer in a wire screen mesh regenerator. Forooghi 
et al. [15] numerically investigated the steady and pulsatile flow 
and heat transfer in a channel lined with two porous layers subject 
to constant wall heat flux under local thermal non-equilibrium 
conditions. Pathak [16] carried out a systematic experimental and 
CFD-based study for the quantification of Darcy permeability and 
Forchheimer coefficients for porous structures under steady and 
periodic flow conditions. 

The aim of this study is to present a stepwise numerical model¬ 
ling approach to obtain the inputs required (viscous/inertia resis¬ 
tance and heat transfer coefficient) for accurate modelling of the 
regenerator as a porous media (Fig. 1). 


These first two steps [17,18] involve the validation of the 
numerical pressure drop and heat transfer characterization, 
through a small detailed 3-D stacked woven wire regenerator 
matrix with the use of well-known experimental correlations 
obtained under oscillating flow conditions [19 . Considering that 
the numerical results agree well with the empirical correlations, 
the approach is further extended to wound woven wire matrices 
to numerically derive the friction factor and the Nusselt correlation 
equations at local instantaneous Reynolds number [17,18 . 

Finally, the characterization is extended to a simplified equiva¬ 
lent porous media model of a representative detailed 3-D matrix in 
order to model a full length regenerator matrix under oscillating 
regeneration cycles based on local instantaneous friction factor 
and Nusselt number correlation equations. 


2. Numerical modelling study 


2.2. Previous numerical studies: validation 

The main physical phenomena to be considered in the charac¬ 
terization of the Stirling regenerator are the pressure drop and 
the heat transfer and the trade-off between both through the 
regenerator is crucial. 

The first step of the numerical characterization of Stirling 
regenerator is the pressure drop characterization study [17 . A 
detailed 3-D isothermal model with an initial flow field is proposed 
for determining pressure drop characteristics (Fig. 2a), thus provid¬ 
ing a straightforward insight in how to determine frictional losses 
at different geometric configurations in the typical range of Rey¬ 
nolds number for Stirling regenerators (Fig. 2b). The friction factor 
is calculated at local instantaneous Reynolds number. 

The heat transfer characterization study [18] is the second step 
of the numerical characterization of woven wire regenerator under 
non-isothermal flow computations concerning the heat transfer 
effects of pressure drop phenomena. The model complexity is 
increased because the flow simulation is transient (time-depen- 
dent) and non-isothermal due to inclusion of heat transfer effects 
(Fig. 3a). 




Fig. 1 . Methodological approach for Stirling regenerator study. 
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Fig. 2. Pressure drop numerical validation [22]: (a) stacked woven wire matrices: 
friction factor vs. Reynolds number and (b) velocity vector through stacked woven 
wire matrix SI 10 for Re « 50 (0.9 m/s). 


2.2. Stirling engine regenerator modelled as porous media 

The flow through a porous media is modelled by adding an 
extra source term to the momentum equation and another one, 
which accounts for heat transfer in the interface, to the energy 
equation depending on which thermal porous media model is 
employed. In general, the Darcy’s law is used in porous media to 
describe slow or creeping flows and the Forchheimer’s law is used 
for the description of high velocity flows. 

The Darcy’s law is obtained from performing steady-state unidi¬ 
rectional flow experiments for a uniform sand column and can be 
written as Eq. (1): 

Vp = - —v (1) 

a 

The Forchheimer’s law is based on Darcy’s law including the 
inertial effects. The inertial effects start to dominate the flow 
through porous media in the high velocity regime. The 
Forchheimer equation is given as follows: 

v P = -^v-/?p / |v|v i (2) 

The Darcy’s law and the Forchheimer’s law coefficients are both 
commonly obtained experimentally in order to describe the flow 
through a porous medium. Here both coefficients are obtained 
based on previous validated numerical results through a small 
detailed 3-D woven wire matrix domain as a representative of 
the whole Stirling regenerator [17]. The micro scale model includes 
individual wires trying to catch any non-uniformity of 3-D flow 
inside the matrix as it is shown in Fig. 2a. 



(b) 


25 


3 


20 


15 


10 


0 


Nusselt Gedeon 64% 
Nusselt Gedeon 68% 

° All Misaligned Stacked 


e 

° 6' 
Or 

.6 

„ / °/ 
o O / 

I 

/ 

° / 

/ 

o / 

o / 
oo / 


08 


OO o' 

O 00 8 / 


e 


° eo ,'o 
o e 


0 6 # 




8>° 

•.-■ ■j--" ' *” 


°%T 

o cP 

° ° o°<% 0 

o & 8 ®—' 

o<8 o " X 


4 


40 


400 


Re 


Fig. 3. Heat transfer numerical validation [23 : (a) stacked woven wire matrices: 
Nusselt vs. Reynolds number and (b) fluid temperature through stacked woven wire 
matrix SI 10 for Re « 50 (0.9 m/s) at 0.02 s. 


Ergun [20] proposed an expression for the Forchheimer coeffi¬ 
cient based on the experimental investigation of fluid flow through 
packed columns and fluidized beds: 





From Eqs. (2) and (3), the Ergun equation is obtained: 



Vp = -^-^- |v|Vi (4) 

The Ergun coefficient is strongly dependent on the flow regime. 
For the slow flows, the inertial effects are very small and can be 
neglected. This reduces the Forchheimer equation to the Darcy 
equation. As the flow velocity increases, inertial effects also 
increase and the flow adapts to the Forchheimer flow regime. 
These inertial effects are accounted for by the Ergun coefficient 
and the kinetic energy of the fluid. 

Ergun [20], also determined the general form friction factor 
equation and this form is used for many friction factor correlation 
equations in Stirling regenerators [17,21-23]: 



fli 

Re 


+ a 2 



Considering the definition of the pressure drop as follows: 


AP = C/ ^ 2 U “ (6) 

The maximum flow velocity i/ max is obtained by dividing the 
frontal maximum velocity by matrix volumetric porosity. The 
hydraulic diameter of the regenerator matrix, d h} is defined by 
Gedeon and Wood [19] for a woven wire matrix as: 
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d w • n v 
= (i -n„) 

Therefore, Reynolds number can be given as: 



Re 



tfinax 



Considering one-dimension, Eq. (6) is transformed into Eq. (9) 
using Eqs. (5) and (8) where the coefficients of Eq. (4) can be easily 
obtained from two-parameters friction factor correlation equations 
(Eq. (5)). 



/LL ■ CL\ 


L 

~2 ^max 
d h 


+ 


a 2 l 2 

2 d h PfUmax 



The above two parameters based friction factor correlation 
equation (Eq. (5)) is derived here to obtain the source term 
included into momentum equations for two specific wire woven 
wire matrix configurations based on the numerical results. 

In the heat transfer model through a porous media, there are 
two possible thermodynamic scenarios: local thermal equilibrium 
and local thermal non-equilibrium. In the present study, the sce¬ 
nario is the local thermal non-equilibrium, and the heat transfer 
coefficient, which is required to define the source term in the 
energy equations, is obtained based on numerical results. 


2.3. Numerical pressure drop characterization 

In this section the detailed numerical results for two specific 
matrix configurations with a volumetric porosity of 63% and 
110 pm wire diameter, S110-63% (stacked) and W110-63% 
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Fig. 5. Pressure drop (Pa) vs. inlet velocity (m/s) trend-line for volumetric porosity 
of 63% at constant working fluid properties. 


(wound) are presented. These results are used as the input of the 
porous media model for the Stirling regenerator. 

Fig. 4 shows the pressure drop results as a function of the length 
through the regenerator matrix, for stacked and wound configura¬ 
tion for Re « 65 (inflow velocity of 1 m/s). It is observed that the 
pressure drop behaviour is stepwise and uniform through the 
regenerator matrix and the trend is linear. The lower pressure drop 
gradient is obtained through the stacked matrix configuration. 

Fig. 5 presents the plot of the pressure drop as a function of the 
inlet velocity through a stacked and wound matrix configuration. It 
is observed that the trend-line is approximated to a second order 
polynomial which corresponds to Forchheimer’s law (Eq. (2)). 
These numerical results can be extrapolated to determine the coef¬ 
ficients for the equivalent porous media. It should be noted here 
that all numerical results are obtained under constant working 
fluid properties of density and viscosity. 

Fig. 6 shows the relationship between friction factor and 
Reynolds number for the numerical results obtained for a stacked 
and wound matrix configuration. Based on the numerical results, 
a specific friction factor correlation equation is obtained for each 
matrix configuration. The specific correlation equation for 
SI 10-63% is given by Eq. (10): 

1118 

C/ = -^-+l-85 (10) 

And the specific friction factor correlation equation for 
W110-63% is given by Eq. (11): 
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Fig. 4. Pressure drop (Pa) vs. Position (mm) through woven wire matrix of 63% 
volumetric porosity for Re « 65 (1 m/s): (a) stacked and (b) wound. 


Fig. 6. Specific friction factor vs. Reynolds number for stacked and wound woven 
wire matrices of 63% volumetric porosity. 
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Table 1 

SI 10-63% and W110-63% regenerator matrix parameters. 
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Matrix 

n v (%) 

<P (1/flm) 

d h (pm) 

L (mm) 

1/a 

c 2 

SI10-63% 

63 

0.013 

187 

1.41 

2.47E + 9 

2.49E + 4 

W110-63% 

63 

0.013 

187 

1.63 

3.72E + 9 

2.69E + 4 



165.5 
~Re~ 


+ 2.04 



Although the previously proposed three-parameters based gen¬ 
eral correlation equations by Costa et al. [17] can be used for 
stacked and wound woven wire matrix configurations, the above 
derived specific correlation equations (Eqs. (10) and (11)) can be 
directly used for each matrix configuration to model the porous 
media. In both configurations, the two-parameter Ergun form 
(Eq. (5)) fits very well with the numerical results and the introduc¬ 
tion of the third parameter is not required as shown in Fig. 6. 

Considering the above specific friction factor correlation equa¬ 
tions (Eq. (10) and (11)), by rearranging Eq. (9) and making use 
of the relation between maximum velocity and inlet velocity, the 
viscous and the inertial resistance coefficients can be also obtained. 
Table 1 summarizes the resistance coefficients and the matrix 
characteristics for both configurations. 


2.4. Numerical heat transfer characterization 

This section presents the numerical heat transfer results 
obtained from the second step solutions of two specific matrix con¬ 
figurations, SI 10-63% and W110-63%. These results are used as 
the input of the porous media model for the Stirling regenerator. 
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Fig. 8. Wire and working fluid temperature (K) vs. Position (mm) through woven 
wire matrix of 63% volumetric porosity for Re « 65 (1 m/s) at 0.02 s: (a) stacked and 
(b) wound. 


Fig. 7 shows the wire temperature distribution as a function of 
the length through the regenerator matrix for a stacked and wound 
woven wire matrix. A clear stepwise temperature gradient is 
observed in the SI 10-63% matrix while a smooth temperature gra¬ 
dient mainly due to the longitudinal conduction through the wires 
is observed in the W110-63% matrix. However, in both cases the 
temperature across the regenerator matrix is not linear, and the 
trend-line is logarithmic. 

Fig. 8 shows the wire and working fluid temperature distribu¬ 
tion as a function of the length through the regenerator matrix 
for the stacked and wound woven wire matrices. The figure also 
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Fig. 10. Computational domain for equivalent porous media model. 


plots the trend-line for the wire temperature and for the mass- 
average working fluid, it can be appreciated that the temperature 
difference is rather logarithmic than linear through the matrix 
under the flow condition studied. 

The scatter observed is due to the internal structure of the 
regenerator matrix where can be clearly identify velocity gradient 
regions and low velocity (wake) regions where the heat transfer 
mechanism is negatively affected [17,18 . 

Fig. 9 shows the relationship between Nusselt number and Rey¬ 
nolds number for the numerical results obtained for a stacked and 
wound matrix configuration. Based on the numerical results 
obtained for different solid and working fluid properties, a specific 
Nusselt correlation equation is obtained for each matrix configura¬ 
tion. The specific Nusselt correlation equation for SI 10-63% is 
given by Eq. (12): 


Nu = 1.91 +0.17Re 080 (12) 

And the specific Nusselt number correlation equation for 
W110-63% is given by Eq. (13): 


Nu = 2.15 + 0.07Re° 88 (13) 

The Nusselt number correlation equations are used to define the 
heat transfer coefficient through a local thermal non-equilibrium 
porous media model for the full Stirling regenerator. The heat 
transfer coefficient is defined as: 


h 


Nu ■ kf 

d h 



Although these coefficients can be experimentally obtained it is 
known that the experimental studies are limited and expensive 
due to the internal structures of the regenerator matrix. Therefore, 
the numerical approach presented is a cost effective tool compared 
to experimentation that can be used with confidence in the 
characterization of woven wire regenerator matrices to reduce 
experimental cost and improve internal geometry. 
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Fig. 11. Equivalent porous media pressure drop (Pa) vs. inlet velocity (m/s) trend- 
line for volumetric porosity of 63% at constant working fluid properties. 


3. Computational principles 

3A. Governing equations for porous media 

3AA. Momentum equations for porous media 

The porous media is modelled by the addition of a source term 
to the momentum equation. The source term is composed of two 
parts: a viscous loss term (Darcy, the first term on the right-hand 
side of Eq. (15)), and an inertial loss term (the second term on 
the right-hand side of Eq. (15)). 

Si = - ( Ypnm + f^ApfWVj) (is) 

\j =i j =i / 

This momentum sink contributes to the pressure gradient in the 
porous cell, creating a pressure drop that is proportional to the 
fluid velocity (or velocity squared) in the cell. Eq. (16) recovers 
the case of simple homogeneous porous media: 

Si = ~{a Vi + C 2 l p f |v|Vi ) (16) 

Assuming an isotropic porous media and based on physical 
velocity, the momentum conservation equation is as follows: 
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Fig. 13. Temperature contours (K) in a section plane a 1 mm for 3-D detailed matrix of 63% volumetric porosity for Re « 65 (1 m/s) at 0.02 s: (a) stacked and (b) wound. 



+ V • (yPfVV ) = -yVp + V • (y f) + yB f 




3A 2. Energy equations for porous media 

For the local thermal non-equilibrium porous media approach, 
a solid zone that is spatially coincident with the porous fluid zone 
is defined, and this solid zone only interacts with the fluid with 
regard to heat transfer. The conservation equations for energy 
are solved separately for the fluid and solid zones. The conserva¬ 
tion equation solved for the fluid zone is: 


d 

Ft 




d 

dXi 


^di +TijVj 


Z h Ji 


+ S h f + h fs A fs (T f - T s ) 


and the conservation equation solved for the solid zone is: 


(18) 


d d 

Q t ^-y)p^=^ 


(\-y)k s ^j+S h s -h fs A fs (T f -T s ) (19) 


32. Numerical methodology 

The fluid is considered to be viscous, incompressible, Newtonian 
and 2-D through the equivalent porous media zone using true 
(physical) velocity conditions under laminar flow assumption. The 
3-D detailed geometry is relevant for the modelling of a woven wire 
regenerator matrix because the flow characteristic is intrinsically 
3-D. However, modelling an equivalent isotropic porous media a 
simplified 2-D domain fully represents the detailed 3-D model. 
On the other hand, the 2-D domain is the optimal model for 
minimizing computational memory and time which is especially 
important under oscillating flow conditions. 

The governing flow equations are discretized and solved 
sequentially using a finite volume method (FVM) based numerical 
flow solver [24 with a second order upwinding scheme for the dis¬ 
cretization of the continuity, momentum and energy equations. 
The convergence criterion for all the velocity components and for 
the continuity is set to 10 s , whereas for energy it is set to 10 8 
for all simulations. The non-dimensional time step, t = u max t/s is 
expressed as the product of the maximum velocity component at 
the inflow boundary and time elapsed divided by the volume cell 
size, in the porous media study the non-dimensional time step val¬ 
ues are chosen as 1. 
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Fig. 14. Pressure drop (Pa) vs. crank angle (°) through: (a) stacked and (b) wound 
matrices. 



32.1. Numerical methodology for equivalent porous media 

The flow is unidirectional, unsteady and all fluid properties 
(including density, viscosity, specific heat, and conductivity) are 
assumed constant due to the small range of temperature used in 
each case (micro-scale model). 


Fig. 15. Temperature profile (K) vs. Position (mm) through porous matrix for: (a) 
stacked and (b) wound. 

The present simplified porous media model simulations are 
carried out by considering the following simplified boundary 
conditions: 


322. Numerical methodology for full regenerator as a porous media 

For the simplified full regenerator solution, the flow is oscillat¬ 
ing, unsteady, and the properties of the working fluid and the solid 
metal wire depend on the temperature. The working fluid is 
nitrogen and the metal wire matrix material is stainless steel. 

3.3. Computational domain and boundary conditions 

3.3.1. Computational domain and boundary conditions for equivalent 
porous media 

Fig. 10 illustrates the 2-D domain of the flow of interest (geom¬ 
etry set-up) which includes a porous zone. The porous zone is 
equivalent to the detailed woven wire matrix previously character¬ 
ized by the authors [17,18]. The rest of the domain remains the 
same and the boundary conditions and simulation options are 
not modified with respect to previous heat transfer studies for 
woven wire matrix, SI 10-63% and W110-63%. 

The effect of the mesh resolution and time step size on the pres¬ 
ent flow is assessed through a grid independence study for three 
different mesh systems containing non-uniformly hexahedral grid 
cells and for three non-dimensional time step (from 0.1 to 1) at a 
Reynolds number of 200. In the temperature distribution values 
obtained from different mesh configurations and a non- 
dimensional time step, there is no major difference among the 
computed values through a porous media. 


1. Inflow boundary: Simplified uniform velocity inlet boundary 
conditions at constant temperature are used to define the fluid 
uniform velocity profile. 

2. Outflow boundary: Pressure outlet boundary conditions are 
assigned to define the static (gauge) pressure by eliminating 
the reverse flow problem. 

3. Side boundary: Free-slip symmetry flow boundary conditions at 
the two side boundaries of the computational domain are 
imposed. 

4. Porous media zone: The settings and inputs for the porous media 
are the following: physical velocity formulation, viscous and the 
inertial resistance coefficients summarized in Table 1, isotropic 
porosity of 63%, interfacial area density (see Table 1) and con¬ 
stant heat transfer coefficient. 

3.3.2. Computational domain and boundary conditions for full 
regenerator as a porous media 

The idea is to model a simplified full size regenerator porous 
media under oscillating regeneration cycles. The full size regener¬ 
ator matrix consists of the total length of a Stirling regenerator 
used in a micro-CFlP reference engine (30 mm), and the full flow 
area in the regenerator is not modelled. The computational domain 
is constructed of hexahedral meshes and 200 time steps per cycle 
at 25 FIz is chosen. 

The porous media model simulations are carried out by consid¬ 
ering the following boundary conditions: 
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Fig. 16. Temperature fluctuations (K) in three regenerator sections vs. time (s) for: 
(a) stacked and (b) wound. 



Fig. 17. Temperature fluctuations (K) in the mid-section vs. time (s) for: (a) stacked 
and (b) wound. 


1. Inflow boundary: Oscillating inlet velocity boundary conditions 
at constant temperature (873 K, corresponding with the engine 
hot end) is used to define the fluid sinusoidal velocity profile 
defined with a user-defined-function (UDF) as follows: 

u = ii 0 sin(27t/t) (20) 

where u 0 is the amplitude of the oscillating inlet velocity (1.5 m/s) 

and/is the engine frequency (25 Hz). 

2. Outflow boundary: Pressure outlet boundary conditions at con¬ 
stant temperature (423 K, corresponding with the engine cold 
end). The working pressure is constant at 26 bar. 

3. Side boundary: Free-slip symmetry flow boundary conditions at 
the two side boundaries of the computational domain are 
imposed. The normal velocity components and the normal gra¬ 
dients of all velocity components are assumed to have a zero 
value. 

4. Porous media zone: The settings and inputs for the porous media 
are the following: physical velocity formulation, viscous and the 
inertial resistance coefficients summarized in Table 1, isotropic 
porosity of 63%, constant interfacial area density (see Table 1) 
and variable heat transfer coefficient. The variable heat transfer 
coefficient of the porous medium is defined through a UDF 
using Eqs. (8), (12), (13) and (14). 

4. Results and discussion 

4.1. Numerical validation for equivalent porous media 

The results obtained from the numerical simulations, which are 

performed for the equivalent 2-D porous media of the stacked 


woven wire matrix SI 10-63% and of the wound woven wire matrix 
W110-63% are described for the detailed matrix configuration. 

Fig. 11 shows the pressure drop results as a function of the inlet 
velocity for SI 10-63% and W110-63% equivalent porous media. 
The figure also shows the numerical results obtained from the 
detailed regenerator matrix simulation solutions for both configura¬ 
tions. Excellent agreement between the results of 3-D regenerator 
matrix and the equivalent 2-D porous media are obtained using 
the viscous and inertia resistance coefficients derived from the 
detailed 3-D regenerator matrix summarized in fable 1. 

Fig. 12 shows the numerical results for fluid and solid temper¬ 
ature distribution as a function of the length through the equiva¬ 
lent porous media. Fig. 12a compares the results for the detailed 
stacked woven wire matrix against its equivalent fluid and solid 
porous media, and Fig. 12b for the wound woven wire matrix. 
The fluid temperature distribution for the detailed matrix corre¬ 
sponds to the mass-average temperature. It is observed from the 
figure that there is a good agreement between the detailed matrix 
temperatures and the equivalent porous media. It is noteworthy 
here that the temperature distribution through the isotropic por¬ 
ous media is smooth and the present model cannot reproduce 
the stepped behaviour in the solid wire matrix solution. 

Fig. 13 represents the temperature contours in a cutting plane at 
1 mm through both detailed 3-D regenerator matrices. The figure 
signifies that the temperature gradient in the fluid is significant 
for both configurations as previously observed in Fig. 6. Through 
an isotropic porous media in a transversal cutting plane, the tem¬ 
perature of the fluid and solid are constant. For this reason the dif¬ 
ferences observed in the results of the equivalent porous media and 
the 3-D model are acceptable considering the temperature gradient 
behaviour observed in the figure. Therefore, it is believed that the 
results verify the suitability of the local thermal non-equilibrium 
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Fig. 18. Porous matrix friction factor and specific friction factor correlation 
equation vs. crank angle (°) for: (a) stacked and (b) wound matrices. 

porous media model as an equivalent of the detailed 3-D regenera¬ 
tor matrices for pressure drop and heat transfer phenomena. 

4.2. Numerical validation for full regenerator as a porous media 

In this section, the numerical simulations for the simplified full 
length 2-D equivalent porous media of the stacked and of the 
wound woven wire matrix are performed under oscillating flow 
conditions. In these models under oscillating regeneration cycles 
the desired solution is the quasi-steady solution, when the results 
for two consecutive cycles are identical. The results present and 
discussed here are obtained in a quasi-steady state for the working 
pressure of 26 bar and the flow oscillating frequency of 25 Hz. 

Fig. 14 shows the variations of the computed pressure drop 
through the porous media matrix for both configurations during a 
complete quasi-steady cycle. It is noted that for the frequency stud¬ 
ied, the pressure drop varies sinusoidal and a small phase lag during 
the first half of the cycle (hot-blow) with respect to the inlet fluid 
velocity is observed. The cycle-averaged pressure drop for the 
SI 10-63% equivalent regenerator porous media is 4700 Pa and for 
the W110-63% equivalent regenerator porous media is 5780 Pa. 

The temperature profiles through the regenerator porous 
matrix for the fluid and solid are presented in Fig. 15. A logarithmic 
and linear trend-line are both included into the figure as reference 
lines. It can be found that the temperature profile through the 
regenerator can be better fitted to a logarithmic curve (regression 
coefficient of about 1). rather than linear line. Therefore, the 
logarithmic mean temperature difference between fluid and solid 
is used for the present heat transfer calculations. 

Fig. 16 illustrates the temperature fluctuations in three sections 
of the regenerator porous media matrix (mid-plane, hot-end 
(1 mm) and cold-end (1 mm)) during three complete quasi-steady 




Fig. 19. Porous matrix Nusselt number and Specific Nusselt number correlation 
equation vs. crank angle (°) for: (a) stacked and (b) wound matrices. 

cycles. It is observed that in different axial sections the solid matrix 
temperature corresponds to the mean fluid matrix temperature. 
On the other hand, both fluid and solid temperatures change peri¬ 
odically with the oscillating flow and the variations of the fluid 
temperature is larger than the solid temperature (almost constant). 

The temperature fluctuations in the mid-section for stacked and 
wound matrix configurations are also illustrated in Fig. 17. A phase 
lag is clearly identified for the matrix temperature due to the heat 
transfer process. Basically, the fluid temperature is higher during 
the first half of the oscillating cycle, while at the second half, the 
solid temperature is higher. 

Fig. 18 represents a comparison between the friction factor 
numerically obtained from the 3-D detailed matrix (Eqs. (10) and 
(11)) at local instantaneous Re number with the friction factor 
obtained through a whole porous media under oscillating flow 
conditions at average Re number. The figure clearly shows that 
the profile of the computed local instantaneous friction factor cor¬ 
responds well with that of average friction factor. 

Fig. 19 illustrates the Nusselt number variation with respect to 
the crank angle for the specific numerically obtained correlation 
from the 3-D detailed matrix (Eqs. (12) and (13)) with the average 
Nusselt number obtained in the present study through a porous 
media under oscillating flow conditions. The Nusselt number peaks 
are due to the very small temperatures difference between the 
fluid and the solid matrix during direction change under oscillating 
flow. The mean Nusselt number obtained through the porous 
media under oscillating flow conditions characterize the heat 
transfer performance of the 3-D detailed regenerator matrix. 
Comparing the temperatures results shown in Fig. 16, there is 
not a significant variation in the temperatures obtained for both 
configurations. However, if the thermal efficiency is calculated a 
noticeable difference may be obtained. 



























































































































S.C. Costa et all Energy Conversion and Management 89 (2015) 473-483 


483 


5. Conclusions 

The present study proposes a stepwise numerical procedure for 
the derivation of the flow resistance coefficients and thermal non¬ 
equilibrium heat transfer coefficient for porous media of Stirling 
regenerator under unidirectional flow conditions. The accuracy of 
the equivalent porous media modelling technique is verified 
against the validated numerical results obtained for the detailed 
3-D woven wire matrices. The methodology is successfully adopted 
for the flow and heat transfer characterization of two different 
woven wire matrix configurations, stacked and wound with 63% 
of volumetric porosity. 

It is remarkable that for specific matrix configuration, the fric¬ 
tion factor correlation equation form that better fits with the 
numerical results is a two-parameter Ergun form. In the case of a 
general correlation equation a third parameter could be necessary 
to better fit the curve at high Reynolds numbers [17,19]. On the 
other hand, the heat transfer results demonstrate that the trend- 
line obtained for the temperature profiles through the regenerator 
porous matrix for the fluid and solid can be better fitted to a loga¬ 
rithmic curve rather than linear under oscillating flow conditions. 

The present results suggest that the pressure drop and heat 
transfer performance of the different configuration matrices can 
be evaluated and compared as well as the thermal efficiency or 
the figure of merit. In the characterization of a stacked and a 
wound woven wire matrices is observed a better overall perfor¬ 
mance in terms of pressure drop and heat transfer for the stacked 
configuration. 

It is believed that the numerical characterization methodology 
presented is advantageous for the design and optimization of the 
regenerator matrix under diverse working conditions (oscillating 
frequency, working pressure, working gas, matrix material, etc.) 
and useful for the multi-D Stirling engine models. 
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